clear

global xmean x2mean Deltax2 r dr Amat px x dx

 q1=0;
 q2=0;
 q3=0;
 q4=0;
 q5=0;
 q6=0;


dx=1;
xmax=1200;
x=0:dx:xmax;  % pocty fotonu
rmax=sqrt(xmax);
dr=rmax/550;
r=0:dr:rmax;  % amplitudy koherentnich stavu

[R,X]=meshgrid(r,x);
% Conditional probability of X on R
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Amat=exp(-(X-R.^2).^2./(R+.00001).^2/2)/sqrt(2*pi)./(R+.00001); % Gaussian approximation
%aux1=ones(size(Amat,1),1)*sum(Amat)*dx;
%Amat=Amat./aux1;
Amatsiz=size(Amat);
MaxNexact=70; % up to how many photons the probability is calculated as exact Poissonian
[Rp,Xp]=meshgrid(0:dr:sqrt(MaxNexact),0:dx:MaxNexact);
AmatExact=exp(-Rp.^2).*Rp.^(2*Xp)./gamma(Xp+1);
AmEsiz=size(AmatExact);
Amat(1:min(AmEsiz(1),Amatsiz(1)),1:min(AmEsiz(2),Amatsiz(2)))=AmatExact;
%
%figure(1)
%mesh(r,x,Amat.*(Amat<1.1))



% Definice funkci
%%%%%%%%%%%%%%%%%%%%%%
function y=F2(z)
 global xmean x2mean Deltax2 r dr Amat px x dx
 q1=z(1);
 q2=z(2);
 q3=z(3);
 q4=z(4);
% q5=z(5);
% q6=z(6);
 pr=exp(q4*r.^4+q3*r.^3+q2*r.^2+q1*r).*r; 
% pr=exp(q6*r.^6+q5*r.^5+q4*r.^4+q3*r.^3+q2*r.^2+q1*r).*r;
 pr=pr/sum(pr)/dr;
 pxN=(Amat*pr')'*dr;
 xmeanNp=pxN*x'*dx;
 x2meanNp=pxN*(x.^2)'*dx;
 y=sum((pxN-px).^2)*dx; % +(xmeanNp-xmean)^2/Deltax2+0.01*(x2meanNp-x2mean)^2/Deltax2^2;
% y=sum(px.*log(px./(pxN+1e-15)));  % Kullback Leibler
endfunction


%x0=-92.083;
%sigma=178.16;
%x0=-22.283
%sigma=158.21
%x0=8.3513
%sigma=155.85
%x0 =  24.509
%sigma =  164.44
%x0 =  43.691 
%sigma =  175.39 
%x0 = 62.801  
%sigma = 191.54  
x0 =  89.822
sigma =209.67 


px=exp(-(x-x0).^2/2/sigma^2);
norm=sum(px)*dx;
px=px/norm;

sum(px)*dx

xmean2=x*px'*dx
x2mean2=x.^2*px'*dx


figure(1)
plot(x,px)




krit=10000000;
PocetPokusu=3250;
resvals=zeros(1,6);
for n=1:PocetPokusu
 %pomodhad=[-1,1,1,1,-1,-1]
 %pomodhad=[-1,1,1,-1]
% pomodhad=5*randn(1,6);
% pomodhad=0.05*randn(1,4);
% pomodhad=[randn(1,3),-1];
% pomodhad=[-8.4875e-02,   4.9143e-03,   1.5046e-03,  -1.2800e-04]+1e-3*randn(1,4);
% pomodhad=[-8.5263e-02   6.6822e-03   1.1052e-03  -1.0797e-04]+1e-3*randn(1,4);
%pomodhad=[7.6108e-02  -1.1576e-02   2.4128e-04  -1.3608e-05]+1e-3*randn(1,4);
%pomodhad=[2.6393e-02  -4.1276e-03  -1.2663e-04  -8.6584e-06]+1e-3*randn(1,4);
%pomodhad=[-2.3663e-03  -1.0011e-03   2.0870e-05  -2.0653e-05]+1e-3*randn(1,4);  
%pomodhad=[-2.1172e-03  -6.4991e-04   1.9839e-04  -2.4943e-05]+1e-3*randn(1,4); 
%pomodhad=[-3.5010e-03   2.1207e-03  -1.0831e-04  -1.5621e-05]+1e-3*randn(1,4); 
%pomodhad=[-3.3742e-03   8.4463e-04   9.8283e-05  -1.9516e-05]+1e-3*randn(1,4); 
%pomodhad=[-3.5999e-03   3.0878e-03  -1.3889e-04  -1.3079e-05]+1e-3*randn(1,4); 
pomodhad=[-3.9290e-03   3.6250e-03  -1.4821e-04  -1.0487e-05 ]+1e-3*randn(1,4); 

 vals=fminunc('F2',pomodhad);
 F2(vals)
 if F2(vals)<krit
  krit=F2(vals);
  resvals=vals;
  npok=n;
 endif
end 
vals=resvals
krit
npok
%pomodhad=vals;
%vals=fminunc('F2',pomodhad)
%F2(vals)
q1=vals(1);
q2=vals(2);
q3=vals(3);
q4=vals(4);
%q5=vals(5);
%q6=vals(6);
pr=exp(q6*r.^6+q5*r.^5+q4*r.^4+q3*r.^3+q2*r.^2+q1*r).*r;
pr=pr/sum(pr)/dr;


figure(2)
plot(r,pr)

pxN=(Amat*pr')'*dr;

figure(3)
plot(x,pxN)
hold on
plot(x,px,'r')
hold off



%xmeanNp=pxN*x'*dx
%x2meanNp=pxN*(x.^2)'*dx
%[sum((pxN-px).^2)*dx,(xmeanNp-xmean)^2/Deltax2,0.01*(x2meanNp-x2mean)^2/Deltax2^2]

figure(4)
%plot(r,pr./(r+0.0001))
plot(r,q6*r.^6+q5*r.^5+q4*r.^4+q3*r.^3+q2*r.^2+q1*r)
